O’Hara on GitHub


knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.height = 3, fig.width = 12)

library(terra)
library(oharac)
library(tidyverse)
library(here)
source(here('common_fxns.R'))

1 Summary

With impact maps generated for all stressors, means calculated both unweighted across species and FV-weighted across functional entities, we can add the impact layers to generate a comprehensive cumulative human impact map.

2 Data

See scripts 2, 3a, and 3b for vulnerability mapping scripts, and scripts 4a-4d for impact calculation scripts.

3 Methods

3.1 Calculate CHI: spp

Sum impacts across all stressor/vulnerability combos.

imp1_fs <- list.files(here_anx('_output/impact_maps/impact_maps_species'),
                      pattern = 'mean',
                      full.names = TRUE)

imp1_stack <- rast(imp1_fs) %>%
  setNames(str_remove_all(basename(imp1_fs), '.+_spp_|_x_.+|_mean.tif'))
ocean_a_rast <- rast(here('_spatial/ocean_area_mol.tif'))
### Calculate cumulative impacts
chi_spp_raw <- sum(imp1_stack, na.rm = TRUE)
chi_spp <- chi_spp_raw %>%
  mask(ocean_a_rast) %>%
  setNames('chi_spp')

writeRaster(chi_spp, here('_output/cumulative_impact_maps/chi_species.tif'),
            overwrite = TRUE)
### plot with ggplot
chi_spp_plot <- here('figs/chi_species.png')

chi_spp_df <- chi_spp %>%
  as.data.frame(xy = TRUE) %>%
  filter(!is.na(chi_spp))

fill_lbls <- c(0, .1, .25, .5, 1, round(minmax(chi_spp)[2], 2))
xfm <- function(x) {x^.25}
p <- ggplot(chi_spp_df, aes(x, y, fill = xfm(chi_spp))) +
  geom_raster() +
  scale_fill_viridis_c(labels = fill_lbls, breaks = xfm(fill_lbls)) +
  coord_sf() +
  theme_void() +
  labs(title = 'CHI by species mean vulnerability',
       fill  = 'CHI')

ggsave(chi_spp_plot, height = 4, width = 8, dpi = 300)
  
knitr::include_graphics(chi_spp_plot)

3.2 Calculate CHI: functional entity

Sum impacts across all stressor/vulnerability combos.

imp2_fs <- list.files(here_anx('_output/impact_maps/impact_maps_funct_entity'), 
                      pattern = 'mean',
                      full.names = TRUE)

imp2_stack <- rast(imp2_fs) %>%
  setNames(str_remove_all(basename(imp2_fs), '.+_fe_|_x_.+|_mean.tif'))
ocean_a_rast <- rast(here('_spatial/ocean_area_mol.tif'))
### Calculate cumulative impacts
chi_fe_raw <- sum(imp2_stack, na.rm = TRUE)

chi_fe <- chi_fe_raw %>%
  mask(ocean_a_rast) %>%
  setNames('chi_fe')

writeRaster(chi_fe, here('_output/cumulative_impact_maps/chi_funct_entity.tif'),
            overwrite = TRUE)
### plot with ggplot
chi_fe_plot <- here('figs/chi_funct_entity.png')

chi_fe_df <- chi_fe %>%
  as.data.frame(xy = TRUE) %>%
  filter(!is.na(chi_fe))

fill_lbls <- c(0, .1, .25, .5, 1, round(minmax(chi_fe)[2], 2))

p <- ggplot(chi_fe_df, aes(x, y, fill = xfm(chi_fe))) +
  geom_raster() +
  scale_fill_viridis_c(labels = fill_lbls, breaks = xfm(fill_lbls)) +
  coord_sf() +
  theme_void() +
  labs(title = 'CHI by functional entity vulnerability',
       fill = 'CHI')

ggsave(chi_fe_plot, height = 4, width = 8, dpi = 300)

knitr::include_graphics(chi_fe_plot)

3.3 Plot relative difference in CHI

Difference in CHI calculated from FE vulnerability relative to species vulnerability (using only spp included in functional entities): \[\frac{CHI_{FE} - CHI_{spp}}{CHI_{spp}} \times 100\%\]

r <- ((chi_fe - chi_spp) / chi_spp * 100) %>%
  setNames('layer')

r_vals <- values(r); r_vals <- r_vals[!is.na(r_vals)]
r_zlim <- max(abs(quantile(r_vals, .001)), abs(quantile(r_vals, .999)))
values(r)[values(r) > r_zlim] <- r_zlim
values(r)[values(r) < -r_zlim] <- -r_zlim

### set up thresholded rasts as well
r1 <- r2 <- r
r1[chi_fe < .10 & chi_spp < .10] <- NA
r1[r1 > r_zlim] <- r_zlim
r1[r1 < -r_zlim] <- -r_zlim
r2[chi_fe < .25 & chi_spp < .25] <- NA
r2[r2 > r_zlim] <- r_zlim
r2[r2 < -r_zlim] <- -r_zlim
 
### Difference plots
chi_diff_plot  <- here('figs/chi_fe_diff_spp.png')
chi_diff1_plot <- here('figs/chi_fe_diff_spp_0.10.png')
chi_diff2_plot <- here('figs/chi_fe_diff_spp_0.25.png')

chi_diff_df <- as.data.frame(r, xy = TRUE) %>%
  filter(!is.na(layer))

map_cols <- hcl.colors(palette = 'Red-Green', n = 15, rev = TRUE)

p <- ggplot(chi_diff_df, aes(x, y, fill = layer)) +
    geom_raster() +
    scale_fill_gradientn(colors = map_cols, 
                         limits = c(-r_zlim, r_zlim)) +
    coord_sf() +
    theme_void() +
    labs(title = 'CHI difference between FE and species',
         fill = '%∆CHI')
  
ggsave(chi_diff_plot, height = 4, width = 8, dpi = 300)


### Difference cropped by .10 threshold

chi_diff1_df <- as.data.frame(r1, xy = TRUE) %>%
  filter(!is.na(layer)) 
p <- ggplot(chi_diff1_df, aes(x, y, fill = layer)) +
    geom_raster() +
    scale_fill_gradientn(colors = map_cols, 
                         limits = c(-r_zlim, r_zlim)) +
    coord_sf() +
    theme_void() +
    labs(title = 'CHI difference between FE and species (0.10 threshold)',
         fill = '%∆CHI')
  
ggsave(chi_diff1_plot, height = 4, width = 8, dpi = 300)


### Difference cropped by .25 threshold

chi_diff2_df <- as.data.frame(r2, xy = TRUE) %>%
  filter(!is.na(layer))
p <- ggplot(chi_diff2_df, aes(x, y, fill = layer)) +
    geom_raster() +
    scale_fill_gradientn(colors = map_cols, 
                         limits = c(-r_zlim, r_zlim)) +
    coord_sf() +
    theme_void() +
    labs(title = 'CHI difference between FE and species (0.25 threshold)',
         fill = '%∆CHI')
  
ggsave(chi_diff2_plot, height = 4, width = 8, dpi = 300)

knitr::include_graphics(chi_diff_plot)

knitr::include_graphics(chi_diff1_plot)

knitr::include_graphics(chi_diff2_plot)

3.4 Table of mean values per stressor

In and outside of neritic zone.

neritic_rast <- rast(here('_spatial/bathy_mol_neritic.tif'))
mean_str_spp <- sapply(imp1_stack, 
                        FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
  setNames(names(imp1_stack))
mean_str_spp_coastal <- sapply(imp1_stack %>% mask(neritic_rast), 
                                FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
  setNames(names(imp1_stack))

mean_str_fe <- sapply(imp2_stack, 
                        FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
  setNames(names(imp2_stack))
mean_str_fe_coastal <- sapply(imp2_stack %>% mask(neritic_rast), 
                                FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
  setNames(names(imp2_stack))

mean_str_df <- data.frame(stressor = names(imp1_stack)) %>%
  mutate(mean_impact_fe = mean_str_fe[stressor],
         mean_impact_fe_nonzero = mean_str_fe_coastal[stressor],
         mean_impact_spp = mean_str_spp[stressor],
         mean_impact_spp_nonzero = mean_str_spp_coastal[stressor]) %>%
  arrange(desc(mean_impact_fe))
knitr::kable(mean_str_df)
stressor mean_impact_fe mean_impact_fe_nonzero mean_impact_spp mean_impact_spp_nonzero
sst_extremes 0.0338281 0.0488863 0.0351647 0.0491788
ocean_acidification 0.0330115 0.0567399 0.0244016 0.0547503
sst_rise 0.0270853 0.0242258 0.0228888 0.0234100
uv_radiation 0.0152301 0.0306834 0.0172201 0.0303766
direct_human 0.0092048 0.0092048 0.0087021 0.0087021
biomass_removal 0.0078020 0.0241021 0.0096138 0.0273267
sea_level_rise 0.0057294 0.0057294 0.0058910 0.0058910
bycatch 0.0040566 0.0155488 0.0042723 0.0151240
shipping_large 0.0023091 0.0108484 0.0025661 0.0121029
light 0.0010080 0.0100309 0.0010055 0.0099937
nutrient 0.0008922 0.0090952 0.0008932 0.0091172
fishing_benthic_dest 0.0001761 0.0016030 0.0001688 0.0015191
benth_str 0.0000046 0.0000449 0.0000043 0.0000415